! This module defines and instantiates objects
! for a level pool type reservoir. The level
! pool reservoir type inherits input and
! output types from the reservoir base
! module and calls instantiation of these into
! sub-objects. The level pool reservoir type
! also points to types for level pool properties
! and state and calls instantiation of these into
! sub-objects. This module also contains the
! subroutine to run level pool reservoir that is
! derived from the reservoir base type interface
! to run reservoir. Running level pool will
! then call the LEVELPOOL_PHYSICS subroutine, which
! processes the given inputs, properties, and
! state for a particular level pool reservoir and
! returns the output/outflow.

module module_levelpool

    use module_levelpool_properties, only: levelpool_properties_interface
    use module_levelpool_state, only: levelpool_state_interface
    use module_reservoir, only: reservoir, reservoir_input, &
                                     reservoir_output
    use module_hydro_stop, only: HYDRO_stop
    use iso_fortran_env, only: int64

#ifdef RESERVOIR_D
    use module_reservoir_utilities, only: create_levelpool_diagnostic_log_file, &
                                          log_levelpool_diagnostic_data
#endif

    implicit none

    ! Extend/derive level pool type from the abstract base
    ! type for reservoirs.
    type, extends(reservoir) :: levelpool

        ! Define pointers to sub-types / sub-objects to and
        ! held by a level pool reservoir object.
        type (levelpool_properties_interface), pointer :: properties => null()
        type (levelpool_state_interface), pointer :: state => null()

        logical :: pointer_allocation_guard = .false.

    contains

        procedure :: init => levelpool_init
        procedure :: destroy => levelpool_destroy
        procedure :: run => run_levelpool_reservoir

    end type levelpool

contains

    !Level Pool Constructor
    subroutine levelpool_init(this, water_elevation,  &
        lake_area, weir_elevation, weir_coeffecient, &
        weir_length, dam_length, orifice_elevation, orifice_coefficient, &
        orifice_area, max_depth, lake_number, lake_opt)

        implicit none
        class(levelpool), intent(inout) :: this ! object being initialized
        real, intent(inout) :: water_elevation           ! meters AMSL
        real, intent(in)    :: lake_area                 ! area of lake (km^2)
        real, intent(in)    :: weir_elevation            ! bottom of weir elevation (meters AMSL)
        real, intent(in)    :: weir_coeffecient          ! weir coefficient
        real, intent(in)    :: weir_length               ! weir length (meters)
        real, intent(in)    :: dam_length                ! dam length (meters)
        real, intent(in)    :: orifice_elevation         ! orifice elevation (meters AMSL)
        real, intent(in)    :: orifice_coefficient       ! orifice coefficient
        real, intent(in)    :: orifice_area              ! orifice area (meters^2)
        real, intent(in)    :: max_depth                 ! max depth of reservoir before overtop (meters)
        integer(kind=int64), intent(in) :: lake_number   ! lake number
        integer, intent(in) :: lake_opt                  ! bypass lake physics (2 to use pass-through)

        character(len=15)   :: lake_number_string

#ifdef RESERVOIR_D
        ! Create diagnostic log file only for development/debugging purposes
        call create_levelpool_diagnostic_log_file(lake_number)
#endif

        if (this%pointer_allocation_guard .eqv. .false. ) then
            ! try to allocate input
            allocate ( this%input )
            if ( .not. associated(this%input) ) then
                ! if the input structure could not be created, call hydro_stop.
                write(lake_number_string, "(I15)") lake_number
                call hydro_stop("ERROR: Failure to allocate level pool input structure for reservoir " &
                // trim(ADJUSTL(lake_number_string)) // ".")
            else
                ! initialize the input structure
                call this%input%init()
            end if

            ! try to allocate output
            allocate ( this%output )
            if ( .not. associated(this%output) ) then
                ! if the output structure could not be created, call hydro_stop.
                write(lake_number_string, "(I15)") lake_number
                call hydro_stop("ERROR: Failure to allocate level pool output structure for reservoir " &
                // trim(ADJUSTL(lake_number_string)) // ".")
            else
                ! initialize the output structure
                call this%output%init()
            end if

            ! try to allocate properties
            allocate ( this%properties )
            if ( .not. associated(this%properties) ) then
                ! if the properties structure could not be created, call hydro_stop.
                write(lake_number_string, "(I15)") lake_number
                call hydro_stop("ERROR: Failure to allocate levelpool properties structure for reservoir " &
                // trim(ADJUSTL(lake_number_string)) // ".")
            else
                ! initialize levelpool properties
                call this%properties%init( lake_area,  &
                    weir_elevation, weir_coeffecient, weir_length, dam_length, &
                    orifice_elevation, orifice_coefficient, &
                    orifice_area, max_depth, lake_number, lake_opt )
            end if
            this%pointer_allocation_guard = .true.

            ! try to allocate state
            allocate ( this%state )
            if ( .not. associated(this%state) ) then
                ! if the state structure could not be created, call hydro_stop.
                write(lake_number_string, "(I15)") lake_number
                call hydro_stop("ERROR: Failure to allocate state properties structure for reservoir " &
                // trim(ADJUSTL(lake_number_string)) // ".")
            else
                ! initialize levelpool state
                call this%state%init( water_elevation )
            end if
            this%pointer_allocation_guard = .true.
        end if

    end subroutine levelpool_init


    !Level Pool Destructor
    subroutine levelpool_destroy(this)
        implicit none
        class(levelpool), intent(inout) :: this ! object being destroyed
    end subroutine levelpool_destroy


    ! Subroutine for running a level pool reservoir,
    ! which will then call the LEVELPOOL method/subroutine for processing the
    ! inputs and returning the output.
    subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, &
        lateral_inflow, water_elevation, outflow, routing_period, dynamic_reservoir_type, &
        assimilated_value, assimilated_source_file)
        implicit none
        class(levelpool), intent(inout) :: this
        real, intent(in)    :: previous_timestep_inflow ! cubic meters per second (cms)
        real, intent(in)    :: inflow                   ! cubic meters per second (cms)
        real, intent(in)    :: lateral_inflow           ! cubic meters per second (cms)
        real, intent(inout) :: water_elevation          ! meters
        real, intent(out)   :: outflow                  ! cubic meters per second (cms)
        real, intent(in)    :: routing_period           ! seconds
        integer, intent(out):: dynamic_reservoir_type   ! dynamic reservoir type sent to lake out files
        real, intent(out)   :: assimilated_value        ! value assimilated from observation or forecast
        character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value


        ! Update input variables
        this%input%inflow = inflow
        this%input%lateral_inflow = lateral_inflow

        ! Update state variables
        this%state%water_elevation = water_elevation

        call LEVELPOOL_PHYSICS(this%properties%lake_number,                  &
                               this%properties%lake_opt,                   &
                               previous_timestep_inflow,                     &
                               this%input%inflow,                            &
                               this%output%outflow,                          &
                               this%input%lateral_inflow,                    &
                               routing_period,                               &
                               this%state%water_elevation,                   &
                               this%properties%lake_area,                    &
                               this%properties%weir_elevation,               &
                               this%properties%max_depth,                    &
                               this%properties%weir_coeffecient,             &
                               this%properties%weir_length,                  &
                               this%properties%dam_length,                   &
                               this%properties%orifice_elevation,            &
                               this%properties%orifice_coefficient,          &
                               this%properties%orifice_area                  &
        )

        ! Update output variable returned from this subroutine
        outflow = this%output%outflow

        ! Set current inflow to previous_timestep_inflow
        this%input%previous_timestep_inflow = inflow

        ! Update water_elevation variable returned from this subroutine
        water_elevation = this%state%water_elevation

        ! The dynamic reservoir type is always set to 1 for level pool in this module because
        ! it cannot change reservoir types
        dynamic_reservoir_type = 1

        ! The assimilated value would always be a sentinel of -9999.0 because level pool does
        ! not assimilate any external observations or forecasts
        assimilated_value = -9999.0

        ! The assimilated source file would always be an empty string
        assimilated_source_file = ""

#ifdef RESERVOIR_D
        ! Log diagnostic data only for development/debugging purposes
        call log_levelpool_diagnostic_data(this%properties%lake_number, inflow, water_elevation, outflow)
#endif

    end subroutine run_levelpool_reservoir

    ! ------------------------------------------------
    !   SUBROUTINE LEVELPOOL
    ! ------------------------------------------------

    subroutine LEVELPOOL_PHYSICS(ln,lake_opt,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa)

        !! ----------------------------  argument variables
        !! All elevations should be relative to a common base (often belev(k))

        real, intent(INOUT) :: H       ! water elevation height (m)
        real, intent(IN)    :: dt      ! routing period [s]
        real, intent(IN)    :: qi0     ! inflow at previous timestep (cms)
        real, intent(IN)    :: qi1     ! inflow at current timestep (cms)
        real, intent(OUT)   :: qo1     ! outflow at current timestep
        real, intent(IN)    :: ql      ! lateral inflow
        real, intent(IN)    :: ar      ! area of reservoir (km^2)
        real, intent(IN)    :: we      ! bottom of weir elevation
        real, intent(IN)    :: wc      ! weir coeff.
        real, intent(IN)    :: wl      ! weir length (m)
        real, intent(IN)    :: dl      ! dam length(m)
        real, intent(IN)    :: oe      ! orifice elevation
        real, intent(IN)    :: oc      ! orifice coeff.
        real, intent(IN)    :: oa      ! orifice area (m^2)
        real, intent(IN)    :: maxh    ! max depth of reservoir before overtop (m)
        integer(kind=int64), intent(IN) :: ln      ! lake number
        integer, intent(in) :: lake_opt ! reservoir physics options (1: levelpool, 2: passthrough)

        real    :: Htmp                ! Temporary assign of incoming lake el. (m)

        !! ----------------------------  local variables
        real :: sap                    ! local surface area values
        real :: discharge              ! storage discharge m^3/s
        real :: tmp1, tmp2
        real :: dh, dh1, dh2, dh3      ! Depth in weir, and height function for 3 order RK
        real :: It, Itdt_3, Itdt_2_3   ! inflow hydrographs
        real :: maxWeirDepth           !maximum capacity of weir
        !real :: hdiff_vol, qdiff_vol   ! water balance check variables
        !! ----------------------------  subroutine body: from chow, mad mays. pg. 252
        !! -- determine from inflow hydrograph

        Htmp = H   !temporary set of incoming lake water elevation...
        !hdiff_vol = 0.0
        !qdiff_vol = 0.0

        !!DJG IF-block for lake model option  1 - outflow=inflow, 2 - Chow et al level
        !pool, .....
        if (LAKE_OPT == 2) then     ! If-block for simple pass through scheme....
#ifdef RESERVOIR_D
           write(6,*) "LEVELPOOL LAKE_OPT=2, using reservoir passthrough"
#endif
           qo1 = qi1                 ! Set outflow equal to inflow at current time
           H = Htmp                  ! Set new lake water elevation to incoming lake el.

        else if (LAKE_OPT == 1) then   ! If-block for Chow et al level pool scheme

           It = qi0
           Itdt_3   = qi0 + ((qi1 + ql - qi0) * 0.33)
           Itdt_2_3 = qi0 + ((qi1 + ql - qi0) * 0.67)
           maxWeirDepth =  maxh - we

           !assume vertically walled reservoir
           !remove this when moving to a variable head area volume
           sap = ar * 1.0E6

           !-- determine Q(dh) from elevation-discharge relationship
           !-- and dh1
           dh = H - we
           if (dh > maxWeirDepth) then
              dh = maxWeirDepth
           endif

           tmp1 = oc * oa * sqrt(2. * 9.81 * ( H - oe )) !orifice at capacity
           tmp2 = wc * wl * (dh ** (3./2.))  !weir flows at capacity

           !determine the discharge based on current height
           if(H > maxh) then
             discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop
           else if (dh > 0.0 ) then              !! orifice and weir discharge
             discharge = tmp1 + tmp2
           else if ( H > oe ) then     !! only orifice flow
             discharge = oc * oa * sqrt(2. * 9.81 * ( H - oe ) )
           else
             discharge = 0.0   !in the dead pool
           endif

           if (sap > 0) then
              dh1 = ((It - discharge)/sap)*dt
           else
              dh1 = 0.0
           endif

           !-- determine Q(H + dh1/3) from elevation-discharge relationship
           !-- dh2
           dh = (H+dh1/3) - we
           if (dh > maxWeirDepth) then
              dh = maxWeirDepth
           endif

           tmp1 = oc * oa * sqrt(2. * 9.81 * ( (H+dh1/3.) - oe ) )
           tmp2 = wc * wl * (dh ** (3./2.))

           !determine the discharge based on current height
           if(H > maxh) then
             discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop
           else if (dh > 0.0 ) then              !! orifice and weir discharge
             discharge = tmp1 + tmp2
           else if ( (H+dh1/3) > oe ) then     !! only orifice flow,not full
             discharge = oc * oa * sqrt(2. * 9.81 * ( (H+dh1/3.) - oe ) )
           else
             discharge = 0.0
            endif


           if (sap > 0.0) then
              dh2 = ((Itdt_3 - discharge)/sap)*dt
           else
              dh2 = 0.0
           endif

           !-- determine Q(H + 2/3 dh2) from elevation-discharge relationship
           !-- dh3
           dh = (H + (0.667*dh2)) - we
           if (dh > maxWeirDepth) then
              dh = maxWeirDepth
           endif

           tmp1 = oc * oa * sqrt(2. * 9.81 * ( (H+dh2*0.667) - oe ) )
           tmp2 = wc * wl * (dh ** (3./2.))

           !determine the discharge based on current height
           if(H > maxh) then  ! overtop condition, not good!
              discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop
           else if (dh > 0.0 ) then              !! orifice and weir discharge
              discharge = tmp1 + tmp2
           else if ( (H+dh2*0.667) > oe ) then     !! only orifice flow,not full
              discharge = oc * oa * sqrt(2. * 9.81 * ( (H+dh2*0.667) - oe ) )
           else
              discharge = 0.0
           endif

           if (sap > 0.0) then
              dh3 = ((Itdt_2_3 - discharge)/sap)*dt
           else
              dh3 = 0.0
           endif

           !-- determine dh and H
           dh = (dh1/4.) + (0.75*dh3)
           H = H + dh

           !-- compute final discharge
           dh = H - we
           if (dh > maxWeirDepth) then
              dh = maxWeirDepth
           endif

           tmp1 = oc * oa * sqrt(2. * 9.81 * ( H - oe ) )
           tmp2 = wc * wl * (dh ** (3./2.))

           !determine the discharge based on current height
           if(H > maxh) then  ! overtop condition, not good!
              discharge =  tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop
           else if (dh > 0.0 ) then              !! orifice and overtop discharge
              discharge = tmp1 + tmp2
           else if ( H > oe ) then     !! only orifice flow,not full
              discharge = oc * oa * sqrt(2. * 9.81 * ( H - oe ) )
           else
              discharge = 0.0
           endif

           qo1  = discharge  ! return the flow rate from reservoir

        !#ifdef HYDRO_D
        !#ifndef NCEP_WCOSS
        !   ! Water balance check
        !   qdiff_vol = (qi1+ql-qo1)*dt !m3
        !   hdiff_vol = (H-Htmp)*sap    !m3
        !22 format(f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f6.0,2x,f20.1,2x,f20.1)
        !   open (unit=67, &
        !     file='lake_massbalance_out.txt', status='unknown',position='append')
        !   write(67,22) Htmp, H, qi1, ql, qo1, dt, qdiff_vol, hdiff_vol
        !   close(67)
        !#endif
        !#endif

        !23 format('botof H dh orf wr Q',f8.4,2x,f8.4,2x,f8.3,2x,f8.3,2x,f8.2)
        !24 format('ofonl H dh sap Q ',f8.4,2x,f8.4,2x,f8.0,2x,f8.2)


        else   ! ELSE for LAKE_OPT....
         call hydro_stop("Invalid lake option supplied to LEVELPOOL_PHYSICS()")
        endif  ! ENDIF for LAKE_OPT....

        return

    ! ----------------------------------------------------------------
    end subroutine LEVELPOOL_PHYSICS
    ! ----------------------------------------------------------------

end module module_levelpool
